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The bulk of this appendix is the output from a Maple worksheet that contains 
many of the "gritty" details of our second order Schwarzschild calculation. 
A Maple worksheet that reproduces the entire calculation is available for 
download at the GRTensor website. 



oo ■ > restart: 



> grtwO: 

GRTensorll Version 1.70 {R5) 
31 May 1998 

Developed by Peter Musgrave, Denis Pollney and Kayll Lake 
Copyright 1994 ~ 1998 by the authors. 
I Latest version available from : http : // astro. queensu.ca/~ grtensor / 

> mine : 

> read 'myutils.mpl' ; 
^ ' > qload(qschw) : 



Default spacetime = qschw 



, For the qschw spacetime 

' Coordinates 



x{up) 
x^ = [r, 9, 0, t] 
Line element 

{l + Ueh2ir,t)+e^H2ir,t))%l) dr^ 



CiJ). ds= ^ ^ + {ehi{r,t) + e^Hiir,t))%l dr d t 



r 



'X ■ + (1 + ^ (e k(r, t) + £2 K(r, t)) %1) d 



^ ' , 2 • tn\2 /I , 1 



+ smiey (1 + 2 t) + K(r, t)) %1) d ^ 

^ (1 - 2 -) (1 - i (£ /jo(r, t) + £2 H^i^r, t)) %1) dt^ 
r I 

%1 :=3cos(6')2 - 1 

Constraints = [e^ = 0, = 0, = 0, = 0, e'^ ^ 0, = 0, = 0, = 0] 

As can be seen from the form of the metric, the lower case functions are the 
first order perturbations, and the upper case functions are the second order 
perturbations. To begin our analysis we calculate the exact contravariant metric 
tensor, which is then quadratically perturbed using the quadpert() routine. 

> grcalc(g(up,up)) ; 

CPU Time = .394 



21 



> gralter(g(up,up) .quadpert, simplify, factor) ; 
ComponerLt simplification of a GRTensorll object: 

Applying routine quadpert to object g(up,up) 

Applying routine simplify to object g(up,up) 

Applying routine factor to object g(up,up) 

CPU Time = .305 

> grdisplay(g(up,up)) ; 

For the qschw spacetime : 
Contravariant metric tensor 
g{up, up) 

g"" = l((-r + 2m)(-4 + 6/i2(r, t) e cosiOf - 2e h2{r, t) - h2{r, tf + &e^h2{r, tf cos{0f 

-9e^h2{r, tf cos{9f - Ge"^ hi{r, tf cos{9f + 96"^ hiir, tf cos{ef + 6 e"^ H2{r, t)cosi9f 
-2e^H2{r, t) + h,{r, tf))/r 

5 * = -ie(3cos(6l)2 - l){'iehi{r, t) h2{r, t) coti{ef - 3 e hi{r, t) ho{r, t) cos{0f - 2 e Hi{r, t) 
+ ehi{r, t) ho{r, t)-shi{r, t) h2{r, t) -2hi{r, t)) 

g^ « = _i(_4 + 6£ k(r, t) cos{ef -2s k(r, t)-9s'^ cos{ef k(r, tf + 6s'^ k(r, tf cos{ef 
+ 6e'^K{r, t) cos{0f - e'^ k{r, tf-2s^Kir, t))/r'^ 

g^ '^ = i(_4 + 6£ k(r, t) cos{0f -2s k(r, t)-9s^ cos{9f k(r, tf+6s^ k(r, tf cos{0f 
+ 6 K(r, t) cos{0f - s^ k(r, tf - 2 K(r, t))/{r'^ (-1 + cos(6i)) (cos(6') + 1)) 

g* * ^ K{4: + 6sho{r, t)cos{0f - 2sho{r, t)-Qs^ho{r, tf cos{0f +9 s'^ cos{0f ho{r, tf 

+ s^ ho{r, tf + 6s^Ho{r, t) cos{0f + 6 s^ hi{r, tf cos{0f - 9 s"^ hi{r, tf cos{0f 
-s^hi{r,tf -2s^Ho{r,t))/{-r + 2m) 

With these quadratically perturbed contra/co- variant components of the metric 
we can now calculate the perturbed Ricci tensor to second order. 

> grcalcalter (R(dn,dn) , 13) ; 

Simplification will be applied during calculation. 

Applying routine 'Apply constraints repeatedly' to object g(dn,dn,pdn) 
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Applying routine 'Apply constraints repeatedly' to object 
Clir(dn,<iii,cln) 

Applying routine 'Apply constraints repeatedly' to object 
Chr(dn,dn,up) 

Applying routine 'Apply constraints repeatedly' to object R(dn,dn) 

CPU Time = 5.811 
> gralter(R(dii,dii) .expand) ; 

Component simplification of a GRTensorll object: 



Applying routine expand to object R(dn,dn) 

CPU Time = .211 

> grdi splay (R (dn, dn) ) ; 

For the qschw spacetime : 
Covariant Ricci 
R{dn, dn) 

R r r = 86717 words. Exceeds grOptionDisplay Limit 
R r e = 28579 words. Exceeds grOptionDisplayLimit 
R r t = 53004 words. Exceeds grOptionDisplayLimit 
R g = 87057 words. Exceeds grOptionDisplayLimit 
R g t = 20438 words. Exceeds grOptionDisplayLimit 
R (p 4> = 72053 words. Exceeds grOptionDisplayLimit 
R t t = 82672 words. Exceeds grOptionDisplayLimit 

For completeness, the definitions of the Lcgcndre fmiction, P2(cos(6')), that 
we will need to calculate the projections listed in Table 1, and the standard 
Zerilli substitutions, that we will need to eliminate Hi, and, K, appear below. 

> p2 := 5/2*(3*cos(theta) ~2-l)/2; 

p2 :=l^cos(0)2-- 
4 4 

> newH(r,t) := 

> (2*r"2-6*iii*r-3*iii"2) / (r-2*in) / (2*r+3*m) *chi(r , t ) +r"2/ (r-2*in) *eta(r 

(2 — 6 TO r — 3 m^) r^ r]{r, t) 

new (r, j .- (r - 2to) (2r + 3to) ^ r - 2m 

> Kdot(r,t) := 6*(r"2+m*r+m"2)/r"2/(2*r+3*m)*chi(r,t)+eta(r,t) ; 

Kdotir,t):=e ^''^T:''f,^f''^ +Vir,t) 
(2 r + 3 to) 
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Following the prescription set out in Table 1, we can now extract the linear 
contributions to the seven non-trivial Ricci components. The resulting PDEs 
are called leqni,,r. 

> al :=coef f (collectC 

> grcomponent(R(dn,dii) , [theta.theta] ) .epsilon) ,epsilon,l) : 

> gl : =coef f (collect (grcoinponent(R(dn,dn) , [phi, phi] ) .epsilon) .epsilon, 1) : 

> agl:=int (expand ( 

> al-gl/sin(theta)"2)*(2*cot(theta)*diff (p2,theta)+2*(2+l)*p2)*sin(theta) 

> ,theta=0. .Pi) ; 

agl := 12/i2(r, t) - 12 ha{r, t) 

> leqn[l]:= {h[2](r,t) = solve(agl,h[2] (r,t))}; 

lequi := {h2{r, t) = ho{r, t)} 

> a2:=subs(leqn[l] ,al) : 

> g2:=subs(leqn[l] ,gl) : 

> gniiap(R(dii,dii) ,subs,leqn[l] , 'x') ; 

Applying routine subs to RCdn.dn) 

> gralter(R(dii,dii) .simplify, expand) ; 

Component simplification of a GRTensorll object: 



Applying routine simplify to object RCdn.dn) 

Applying routine expand to object R(dn,dn) 

CPU Time = 14.021 

> bl :=coeff (collect (grcomponent(R(dii,dn) , [theta.r]) .epsilon) ,epsilon,l) : 

> leqn[2] :=termsinip( 

> collect (l/3*int(bl*diff(p2,theta)*sin(theta) ,theta=0. .Pi) ,lperts)) ; 

leqn2 := -2 — — — + (— hair, t)) - — k r, t)) + — 

[—r + 2mjr or or —r + 2m 

> cl :=coef f (collect (gr component (R(dn,dn) , [theta,t] ) , epsilon) , epsilon, 1) : 

> leqn[3] :=termsimp( 

> collect(l/3*int(cl*diff (p2,theta)*sin(theta) ,theta=0. .Pi) ,lperts)) ; 

,= . MMW _ (-. + 2„)(^..(M)) _ , a ^„ _ ( I ^(^^ ^„ 

> dl : =coef f (collect (grcomponent (R(dn,dn) , [t ,r] ) , epsilon) , epsilon, 1) : 

> leqn[4] : =termsimp(collect (int (dl*p2*sin(theta) ,theta=0 . .Pi) ,lperts)) ; 

leqn, := 3 ^ - (-^ k(., t)) + ^^^^ - ^ f ^) t)) 
r"' at or r r{—r + 2m) 



24 



> el : =coef f (collect (gr component (R(dn,dn) , [r ,r] ) , epsilon) , epsilon, 1) : 

> leqn[5] : =termsimp(collect (int (el*p2*sin(theta) ,theta=0. .Pi) ,lperts)) ; 

- -3 a (I, ^(M)) - (I. Mr. .,) + ^(iMp . r ' * '» 



{-r + 2m)r 2^dr'^ ^' ^dr"^ " -r + 2m 2 {-r + 2my 
^ho{r,t) (|:k(r,i))(3m-2r) m(f/ii(r,i)) 



— r + 2m (— r + 2m)r (— r + 2m)2 

> f 1 : =coef f (collect (gr component (R(dn,dn) , [t ,t] ) , epsilon) , epsilon, 1) : 

> leqn[6] :=termsimp(collect(int(f l*p2*sin(theta) ,theta=0. .Pi) ,lperts)) ; 

le,n, := -3 ^ ^rn) H^ir, t) _ 1 + 2 .o(M)) _ ^^^^ 

(-^ + 2m)(afip/ii(r,t)) 1 . (_^ + 2m) (|: /io(r, t)) 



r 2'9f2 
m(-r + 2m)(|:k(r, t)) (3m - 2r) /ti(r, t)) 



> leqn [7] : =termsimp (collect ( 

> int((a2+g2/sin(theta)"2)*p2*sin(theta) ,theta=0. .Pi) .Iperts)) ; 

/egn7 := 4k(r, i) + 2 /io(r, i) + (— k(r, i)) (-r + 2 m) r ^^^1' - 2 (— /io(r, t)) {-r + 2 m) 

or'' —r + Zm or 

+ 2 (— k(r, i)) (3 m - 2 r) - 2 (- /ii (r, t)) r 

While this completes our formal analysis of the linear problem, these seven 
equations, Zegni. y, must now be combined in some unobvious ways in order 
to facilitate the reduction of the source terms in the second order calculations. 
These details can be found in the "complete" worksheet at the GRTensor web- 
site. 

In order to determine the second order PDEs {qeqni,,y) we simply repeat the 
linear analysis, but now selecting only the components of the Ricci tensor. 
The second order perturbation functions will, necessarily, appear in these PDEs 
in exactly the same manner as the linear perturbation functions appear in the 
leqrii equations. 

> aal : =simplify(coeff (collect ( 

> grcomponent(R(dn,dn) , [theta.theta] ) .epsilon) , epsilon, 2) ,trig) : 

> ggl : =simplify (expand (l/sin(theta) "2*(coef f (collect ( 

> grcomponent(R(dn,dn) , [phi, phi] ), epsilon) .epsilon, 2))) , trig) : 

> zzl:=simplify(expand(ggl-aal)) ; 

zzl := ^ H2{r, t) - ^ Ho{r, t)-^ hi{r, tf + ^ ho{r, tf + y /ii(r, tf cos{9f - y cos{ef ho{r, tf 

+ 9cos(e)^/io(r, tf + ^i?o(r, t)cos{ef - ^i?2(r, t) cos{ef - 9 hi{r, tf cos{0f 
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> zz2 : =iiit (zzl* (2*cot (theta) *dif f (p2 , theta)+2* (2+1) *p2) *sin(theta) , 

> theta=0. .Pi) ; 

zz2 := -12 F2(r, - y ho{r, tf + y hi{r, tf + 12 Ho{r, t) 

> qeqnl:= {H[2](r,t) = solve(zz2,H[2] (r,t))}; 

qeqnl {ff2(r, t) = ho{r, tf + ^ hi{r, tf + Ho{r, t)} 

> qeqn [1] : =qeqnl ; 

qeqn^ := {7?2(r, t) = -^ho{r, tf + ^ hi{r, tf + Ho{r, t)} 

> grmap(R(dn,dn) , subs, qeqnl, 'x') ; 

Applying routine subs to R(dn,dn) 

> aal :=simplify(coef f (collectC 

> gr component (R(dn,dii) , [theta, theta] ) ,epsilon) ,epsilon,2) ,trig) : 

> ggl : =siinplify (expand (1/ sin (theta) ~2*(coef f (collect ( 

> gr component (R(dn,dn) , [phi,phi] ) ,epsilon) ,epsilon,2))) ,trig) : 

> qeqn [2] :=termsimp (collect ( 

> int((ggl+aal)*p2*sin(theta) ,theta=0. .Pi) ,qperts)) : 

> bbl :=simplif y(coef f (collect( 

> grcomponent(R(dn,dn) , [theta, r] ) .epsilon) ,epsilon,2) .trig) : 

> qeqn [3] :=termsimp (collect ( 

> int(bbl*diff (p2,theta)*sin(theta) ,theta=0. .Pi) ,qperts)) : 

> ccl:=simplify(coeff (collect( 

> gr component (R (dn, dn) , [t , theta] ) .epsilon) .epsilon, 2) .trig) : 

> qeqn [4] :=termsimp (collect ( 

> l/3*int(ccl*diff (p2.theta)*sin(theta) .theta=0. .Pi) .qperts)) : 

> ddl : =simplify(coeff (collect ( 

> gr component (R(dn,dn) . [r,t] ) ,epsilon) ,epsilon,2) ,trig) : 

> qeqn[5] :=termsimp(collect(int(ddl*p2*sin(theta) .theta=0. .Pi) .qperts)) : 

> eel : =simplify(coeff (collect ( 

> gr component (R(dn.dn) . [r,r] ) .epsilon) .epsilon,2) ,trig) : 

> qeqn [6] :=termsimp(collect( 

> expand(int (eel*p2*sin(theta) ,theta=0. .Pi)) ,qperts)) : 

> qsubl:= {diff (H[0] (r,t) ,t)=solve(qeqn[4] ,diff (H[0] (r,t) ,t))}; 

f-i t y~i fi 

qsubl := {- Hoir, t) = -{UH,{r, m + 7 (— Hi{r, t)) - 14 (— H,{r, t))mr-7{- K(r, t)) 

+ 3 hoir, t) ho{r, t)) r^ - 3 h^{r, t) h^{r, t)) + 2 k(r, t) k(r, t)) r^)/r^} 



> qeqn[l]; 



{ff2(r, t) = -i ho{r, tf + i h^{r, tf + Ho{r, t)] 
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> undiff (qeqn[4] ) ; 



Hi [r, t) m 




(I H„(r, 0) - (I K(r, ()) + A (| /.„(•■, <)') 



The derivation of the first ZeriUi equation begins with the elimination of Hq, 



using qeqn^, in the R^.^ component of the Ricci tensor. We then make the previ- 
ously defined ZeriUi substitutions, and make various substitutions involving the 
leqrii equations^]. The resulting expression is called etaeqn in our Maple work- 
sheet. To simplify this rather awkward expression, we employ a new routine, 
mcollect(), that collects specified terms that are multiplied together. This rou- 
tine could, for example, collect together terms that have a factor of, ab, but not, 
say, , or 6^ . A more complete description of this routine appears in Appendix 
A. (Readers who are familiar with the computer algebra environment will, no 
doubt, already see the utility of such a procedure.) 

> etal : =qeqn [5] : 

> eta2:=termsimp (collect (subs (qsubl, etal) ,qperts)) : 

> eta3 : =termsimp (collect (simplify ( 

> subs(H[l] (r,t)=newH(r,t) ,diff (K(r,t) ,t)=Kdot(r,t) , 

> diff (K(r,t) ,r,t)=diff (Kdot(r,t) ,r) ,eta2)) ,zfuncs)) : 

> eta4 : =et a (r ,t) -terinsimp( collect (solve (etaS , eta(r , t) ) ,zf uncs) ) : 

> eta5:=termsimp (collect ( 

> subs( {diff (h[0] (r,t) ,r ,r)=solve(leqn[5] ,diff (h[0] (r,t) ,r,r))},eta4) , 

> zf uncs) ) : 

> eta6:=termsimp (collect ( 

> subs( {diff (h[0] (r,t) ,r)=solve(leqn[2] , diff (h[0] (r,t) ,r)) , 

> diff (k(r,t) ,r ,t)=solve(leqn[4] ,diff (k(r,t) ,r,t)) },eta5) ,zfuncs)) : 

> ug:= {dif f (k(r,t) ,r,r)=termsimp(collect(solve( 

> subs (diff (h[0] (r,t) ,r)=solve(leqn[2] ,diff (h[0] (r,t) ,r)) , 

> subs (Isubl , solve (leqn [6] ,diff (h[0] (r,t) ,r,r))- 

> solve(leqn[5] ,diff (h[0] (r,t) ,r,r)))) ,diff (k(r,t) ,r,r)) ,lperts)) 



> eta7:=termsimp(collect(subs(ug,eta6) .zfuncs)) : 

> etaeqn :=eta7: 

In order to produce a simplified final expression, we extract the source terms 
from etaeqn (by setting the second order perturbation functions to zero) and 
then simplify the result by repeated application of mcollect(). 

> esourcel : =mcollect (expand( 

^The actual substitutions are somewhat arbitrary. We choose to eliminate certain pertur- 
bation functions so as to reproduce the results presented in [1]. 



} 
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> eval(subs(eta(r,t)=0,chi(r,t)=0,etaeqii))) , {h[l] (r,t) ,h[0] (r,t)}) : 

> esource2:=incollect(esourcel, {h[l] (r,t) ,diff (h[l] (r,t) ,t)}) : 

> esource3:=mcollect(esource2, {h[l] (r,t) ,diff (k(r,t) ,r)}) : 

> esource4:=mcollect(esource3, {h [0] (r , t) ,dif f (k(r , t) ,t) }) : 

> esourceS :=incollect (esource4, {dif f (h [1] (r , t) , t) ,dif f (k(r ,t) ,t) }) : 

> esource6:=mcollect(esource5, {k(r ,t) ,h[l] (r ,t)}) : 

> esource7:=mcollect(esource6, {k(r,t) ,dif f (k(r,t) ,t)}) : 

> esource8:=nicollect(esource7, {dif f (k(r,t) ,r) ,diff (h[0] (r,t) ,t)}) : 

> esource9:=incollect(esource8, {k(r,t) ,diff (h[0] (r,t) ,t)}) : 

> esourcel0:=mcollect(esource9, {h[0] (r,t) ,diff (h[0] (r,t) ,t)}) : 

> etasource :=kf actor (esource 10, (r-2*in)/7/(2*r+3*m)) : 



It is now a simple matter to put these results together to arrive at the first 
Zerilli equation for the second order calculation: 

> etaf inal :=eta(r ,t) = 

> termsimp (collect (solve (eval (subs (h[0] (r,t)=0,h[l] (r,t)=0,k(r ,t)=0, 

> etaeqn)) ,eta(r,t)) ,zfuncs) ) -etasource ; 

etafinal := v{r, t) = J^rX{r, t)){-r + 2m) _ 1 ^^^^ _ 2^-, ^_ ^^^^^ ^^^^^ 

+ 2k(r, t) (I- ho{r, t)) + r (|- k(r, t)) (|- ho{r, t)) + 2k(r, t) (|^ k(r, t)) + 4 ^'^''^ 



dr 



m k(r, t)) ho{r, t) ^ k(r, t)) h^{r, t)) 




—r + 2m —r + 2m 

This agrees with the result found in [1] . 

The development of the second order Zerilli "wave equation" follows a similar 
procedure to that of the etafinal result. We first take the time derivative of 
qeqn2, eliminate Hq and Hq using qequi, make the Zerilli substitutions, and 
substitute for ry from our first Zerilli result. We then repeat this process for 
qenqs and add the resulting equations such that the x' term vanishes. The 
result of these operations, which is called chieqn in our worksheet, is essentially 
the raw form of the result we are seeking. This expression is, however, rather 
large. 

> atmpl:=diff (qeqn[2] ,t) : 

> atmp2 : =termsimp (collect (expand ( 

> subs (qsubl, dif f (qsubl,r) ,atmpl)) , {diff (K(r,t) ,t,t,t)})) : 

> atinp3:=terinsiinp(collect(expand(eval(subs(H[l] (r ,t)=newH(r ,t) , 

> diff (K(r,t) ,t)=Kdot(r,t) ,diff (K(r,t) ,t, 

> r)=diff (Kdot(r,t) ,r) , diff (K(r,t) ,t ,r ,r)=dif f (Kdot (r ,t) ,r,r) ,atmp2))) .zfuncs)) : 
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> atmp4 : =termsimp (collect (expand(eval (subs (etaf inal , atmp3) ) ) , zf iulcs) ) : 

> coef a:=simplif y (coef f (atmp4,dif f (chi(r ,t) ,t,r,t),l)); 

coefa := —r^ 

> btmpl : =termsimp (collect (expand( 

> subs(qsubl,diff (qsubl.r) ,dif f (qeqn[3] ,t))) , {diff (H[l] (r,t) ,t,t) , 

> diff (K(r,t),t,t,t)})): 

> btmp2 : =termsimp (collect (expand (eval ( 

> subs(H[l] (r ,t)=newH(r ,t) ,diff (K(r,t) ,t)=Kdot(r,t) , 

> diff (K(r,t) ,r,t)=diff (Kdot(r,t) ,r) , btmpl))) .zfuncs)) : 

> btinp3:=tennsimp(collect(expand(eval(subs(etafinal,btinp2))) , zfuncs)) : 

> coefb:=simplif y(coeff (btmpS.dif f (chi(r,t) ,r,t,t),l)); 



coefb := 3 ■ 



^2 



— r + 2m 

> ftmpl : =termsimp (collect (expand(atmp4/coefa-btmp3/coefb) , 
>{chi(r,t) ,diff (chi(r,t) ,r) ,diff (chi(r,t) ,t) ,dif f (chi (r ,t) ,t,t) , 

> diff (chi (r,t) ,r,r) , diff (chi (r,t) ,r,t) , diff (chi (r,t) ,r,t,t) , 

> diff (chi(r,t) ,r,r,r)})) : 

> f tmp2 :=termsimp( collect (expand (eval ( 

> subs(k(r,t)=0,h[l] (r,t)=0,h[0] (r,t)=0, ftmpl))) , {chi(r,t), 

> diff (chi (r,t) ,r) , diff (chi (r,t) ,r,r) ,dif f (chi (r ,t) ,t,t) , 

> diff (chi(r,t) ,r,r,r)})) : 

> chil :=termsimp(collect( 

> diff (chi (r,t) , t,t)-solve(ftmpl, diff (chi (r,t) ,t,t)) .zfuncs)) : 

> chieqn: =f tmpl : 

> nops (chieqn) ; 

421 

Although, at this point we can at least verify that it reproduces the correct 
linear Zerilli equation when the source terms are set to zero. 

> termsimp(collect(expand(eval(subs(k(r,t)=0,h[l] (r,t)=0,h[0] (r,t)=0, 

> chieqn/coeff (chieqn, diff (chi (r,t) ,t,t) ,1)))) , {chi(r,t) ,dif f (chi (r ,t) ,r) , 

> diff (chi (r,t) ,r,r) ,diff (chi(r,t) ,t,t) ,diff (chi(r,t) ,r,r,r) })) ; 

l))(-r + 2m)^ (-r + 2/,0(3//r'' + ()r,//^+4r^;,;+4r'j,v(/-. t) 



+ 2- 



iS^xir, t)){-r + 2m)m 
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This is just the hnear ZeriUi wave equation with the quadrapole potential, as 
required. 

In order to make contact with the results presented in [1] , we now work with 
the "renormalized" function ^. 

> zetasub:= {chi(r,t) = 

> zeta(r ,t)+2/7*(r"2/(2*r+3*m)*k(r,t)*dif f (k(r,t) ,t)+k(r ,t) ~2) 

}; 

zetasub \ x{r, t) ^ ({r, t) + — + - k r, t) S 

I 7 2r + 3m 7 I 

> zetal : =eval (subs (zetasub, chil) ) : 

> nops (expand(zetal) ) ; 

161 



This substitution has considerably reduced the number of terms in the expres- 
sion. We now eliminate ho in favor of using: 



> psisub:={h[0] (r,t) = 

> (psi(r,t) - r/3*k(r,t))*3*(2*r+3*m)/r/(r-2*ni) +r*diff (k(r,t) ,r)}; 

{ipir,t)~l-rk{r,t)){2r + 3m) g 

psisub := { hoir, t) = 3 ^ — + r (— k(r, t)) ] 

r (r — 2 m) or 



After making this substitution, and a fairly extensive simplification of the source 
terms using the linear perturbation equations^, we arrive at the unsimplified 
answer zetaeqn. 

> nops (zetaeqn) ; 

154 

Here again we make use of the mcollect() routine to simplify our result. 

> f inall : =mcollect (expand (zetaeqn) , {psi(r,t) ,diff (psi(r,t) ,r,r)}) : 

> final2 : =mcollect (f inall , {psi (r , t) ,dif f (psi (r ,t) , t) }) : 

> f inal3 : =mcollect (f inal2 , {dif f (psi (r , t) , t) ,dif f (psi (r,t) ,r,r,r)}) 

> f inal4:=mcollect(f inal3, {dif f (psi(r ,t) ,r) ,diff (psi(r,t) ,t)}) : 

> finals : =mcollect (f inal4, {dif f (psi (r , t) , t) ,dif f (psi (r , t) , t ,r) }) : 

> f inal6 : =mcollect (f inalB , {psi(r ,t) ,dif f (psi(r ,t) ,r)}) : 

> f inal7 : =mcollect (f inal6 , {dif f (psi (r , t) ,r) ,dif f (psi (r , t) , t ,r ,r) }) 

> finals : =mcollect (final? , {dif f (psi (r , t) , t) ,dif f (psi (r , t) ,r ,r) }) : 

> final9 : =mcollect (finals , {psi(r ,t) ,dif f (psi(r ,t) ,r ,t)}) : 



^While we have omitted these details here, they can be found in the full Maple worksheet 
at the GRTensor website. 
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> f inallO:=mcollect(f inaig, {diff (psi(r,t) ,r) ,diff (psi(r,t) ,r,t)}) : 

> f inalll :=incollect(f inallO, {psi(r,t) ,diff (psi(r,t) ,t,r,r)}) : 

> final 12 :=mcollect(f inalll , {diff (psi (r ,t) ,r) ,dif f (psi(r ,t) ,r ,r)}) : 

> f inall3:=mcollect(f iiiall2,{psi(r,t)}) : 

> f inall4:=mcollect(finall3, {diff (psi (r,t) ,t) }) : 

> f inall5:=mcollect(finall4, {diff (psi (r,t) ,r) }) : 

> f inall6 : =mcollect (f inallB , {diff (psi(r,t) ,r,r)}) : 

> f inall7:=mcollect(finall6, {diff (psi(r,t) ,r,t)}) : 

> nops(f inall7) ; 

42 

We now make one final set of substitutions, /i = (r — 2m), and, A = (2r + 3m), 
and extract the source terms from the second order Zerihi wave equation; 

> zetasource:=subs((r-2*m)=mu, (2*m-r)=-mu, (2*r+3*m)=lambda, 

> kf actor (evaKsubs (zeta(r , t)=0 , -f inall7) ) , 12/7* (r-2*m) '3/ (2*r+3*in) ) ) ; 

12 3 / 1 (7m-3r)(f V;(r, t)){^^i;{r, t)) 

zetasource := —u - — 5 ^^^^ 

7 I 3 II 

(8r2 + 12mr + 7m2)V;(r, t){^^P{r, t)) 



+ 



-4 



- 12 



{2r^ + Ar^m + 9rm? + 6m^)i>{r, t)%l 
(-2 - 9 r'' m - 6 m^ + 2 m^ + 15 r m"* + 15 m^) ■i/'(r, t)'^ 



+ 



(J? A 

(112 + 480 m + 692 m? + 762 m^ + 441 r + 144 m^) ip{r, t) {^i;{r,t)) 

fj? A^ 



^ ir' + mr + m')i§jHr,t))i^^i;{r,t)) ^ M%1^A 
fj? 3 
_^ 1 (12r3 + 36 r2m + 59 rm^ + 90 m^) {-§p^{r, t)f 



- 2 



3 /xr^ 
(32 + 88 m + 296 m? + 510 m^ + 561 r + 270 m^) i){r, i) (^ Vl'^, *)) 



/xr7A2 

1 (-18 + 4 m + 33 rm2 + 48 m3)(^V(?-, i))(^V'(^, t)) 
~ 3 /i2 r4 A 

l {lMr,t)){-^i^{r,t)) mi,{T,t){^i,{r,t)) ^ jm^ - 2r^) {§-^^P{r, t))%l 
3 r2 A fxr^ X 

^ 4 (3r^ + 5mr + 6m^)(|:V>(r,t))%l ^ 1 (|, V,(r, t)) (g^ V,(r, t)) 
3 3 r2 



(r^ + mr + m2)2(^V,(r,t))2 1 (^^ ^(r, t))2 A ' 
r"* /x'^ A 3 /X r2 
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%1 :=|^^(r, t) 
This is the effective source term for the second order Zerilh equation. 
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Abstract 

bb" 
0\ . 

^ This article outlines our derivation of the second order perturbations to a Schwarzschild black 

hole, highlighting our use of, and necessary reliance on, computer algebra. The particular per- 
turbation scenario that is presented here is the case of the linear quadrapole seeding the second 
order quadrapole. This problem amounts to finding the second order Zerilli wave equation, and in 
particular the effective source term due to the linear quadrapole. With one minor exception, our 
QQ ■ calculations confirm the earlier findings of Gleiser, et. al. On route to these results we also 

illustrate that, with the aid of computer algebra, the linear Schwarzschild problem can be solved in 
a very direct manner (i.e., without resorting to the usual function transformations), and it is this 
"direct method" that drives the higher order perturbation analysis. The calculations were performed 
using the GRTensorll computer algebra package, running on the Maple V platform, along with sev- 
ly-^ I eral new Maple routines that we have written specifically for these types of problems. Although we 

i have chosen to consider only the "quadrapole-quadrapole" calculation in this article, the GRTensor 

I environment, with the inclusion of these new routines, would allow this analysis to be repeated 

for a far more general problem. These routines, along with Maple worksheets that reproduce our 
calculations, are publicly available at the GRTensor website: www.astro.queensu.ca/~grtensor . The 
interested reader is invited to download and use them to reproduce our results and experiment. 

cr| 

5-1 . Introduction and Overview 

The use of perturbative methods within general relativity to investigate the stability of black holes 
k>( \ has, at the linear level, a long and successful history. The need for such analysis is obvious. While 

' these objects are very mathematically rich and interesting, we need to know if they are indeed 

I physical objects. If it were found that black holes were unstable to small perturbations we certainly 

would not be justified in expecting them to be the ultimate end state of runaway gravitational 
collapse. Of course, the linear analysis has shown that black holes are in fact stable within the context 
of first order perturbation theory. Regge and Wheeler |Q were the first to successfully perform a 
linear black hole analysis in the late 1950's with their ground-breaking study of the (spherically 
symmetric) Schwarzschild metric. Their results were later clarified and refined by Zerilli |^ in 
1970, and it is this work that has become the standard for the linear analysis of the Schwarzschild 
problem. It is important to note that all of these calculations were carried out using the "intuitive" 
metric perturbation (MP) method, i.e., by analyzing perturbations of the metric tensor itself. This 
procedure will be described in greater detail below, but for now we simply note that this stands 
in stark contrast to the only successful linear analysis of the (rotating, axis symmetric) Kerr black 
hole. It is now not difficult to show by direct calculation that when one linearizes the Einstein 
equations about the Kerr vacuum solution the "regular" angular modes, i.e., spherical harmonics, 
do not decouple. This is nothing less than a disaster for the metric perturbation analysis. Here one 
can no longer uniquely project out the contributions from each of the infinite number of multipoles 
required to construct a general perturbation. 

The successful analysis of the linearized Kerr problem, first performed by Teukolsky pT in the 
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early 1970's, avoids this problem of angular mode coupling entirely by recasting the problem in 
terms of the Neuman-Penrose (NP) formalism. (The NP formalism is an example of the more 
general tetrad approach to general relativity, that is distinguished by the fact that all of the tetrad 
vectors are chosen to be null. See for example, Chandrasekhar [5].) Much of the power that the NP 
approach brings to linear perturbation problems is a direct consequence of its formulation of the 
Einstein equations. Here six of the Einstein equations are linear(!) at the outset. An immediate 
result of this rather remarkable fact is that any quantities that are zero in the background solution 
automatically appear as linearized perturbations. The essence of this approach then, is that one does 
not have to linearize the Einstein equations at all. It should be noted that performing a "useful" 
linearization of the Einstein equations about a given metric, even with the aid of modern computer 
algebra packages, can be a non-trivial task. While the NP and MP methods can in principle be 
related to each other via the perturbed tetrad vectors, which must combine to form the perturbed 
metric functions, the two approaches can appear to have surprisingly little in common. We shall 
return to this point in our concluding remarks. As a final comment on the linearized Kerr analysis 
we note that while the NP scheme does make the problem tractable, the analysis remains rather 
complex, and is, in the words of Chandrasekhar j^, "...prolixious in its complexity.". 

The ultimate goal of this article is to present our calculations of second order perturbations to the 
Schwarzschild metric, with an emphasis on the computer algebra methods, packages, and techniques 
that we employ. To properly set the stage for this end result we must first outline the details of 
the linear perturbation theory and its connection to the second, and higher, order perturbations. 
This material appears in the following section. The next section will then give an example of the 
linear theory at work. Here we will follow through the linear Schwarzschild calculation and see 
that with the aid of modern computer algebra this analysis can be carried out in a thoroughly 
transparent fashion. Here we will not derive the usual Zerilli results, but rather illustrate, and find 
the prescription for, the decoupling of the linear perturbation modes. With this material covered we 
will then move on to illustrate the second order calculations. The main goal of this section will be 
to find the form of the effective potential in the Zerilli wave equation. This is where we will see the 
absolute necessity of employing computer algebraic methods to perform these types of calculations. 
We will also see that with the addition of a few new simplification routines the analysis can be made 
much simpler. The final section will, of course, summarize our results and attempt to view them in 
the broader context of more general problems. The current direction of our ongoing research will 
also be discussed. 

Perturbation Theory: A Short Review 

At the outset it is important that we clearly define what we mean by a perturbation, and its relative 
order. Equally important is that we come to fully understand the linear problem. Indeed, once the 
linear problem has been completely solved, we will see that the second order problem is, in principle, 
"just" a matter of computational volume. In what follows we will limit our discussion to the MP 
method, and give only a brief description of the methodology. (A discussion of some of the subtler 
points of the associated questions of gauge can be found in [1].) We start by writing the metric 
as a known, background, vacuum solution with the addition of higher order terms that mimic an 
expansion in a "small" perturbation parameter, e, 

and substitute this "expansion" into the vacuum Einstein equations, 

Rap = . 

To produce an "n*'"' order perturbation we simply truncate all terms with powers higher than 0(e") 
in the resulting equations; so that a first order perturbation involves keeping only up to 0(e) terms. 
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This, along with our restriction to vacuum background solutions, allows us to express the linearly 
perturbed Einstein equations as a homogeneous linear operator, that depends on the background 
metric, acting on the linear perturbations^]: 

e/:^<o, (gW) = . 

The next step in the MP process is to perform the "usual" mode expansion of the perturbation 
functions, 

oo 
i=0 

where g^^^ represents all of the perturbed metric functions, and the Ui{9,(j)) are typically spherical 
harmonics or Legendre functions. If all goes well, when these expansions are put through the linear 
operator the angular modes will decouple. The mathematical statement of this is, 




where /i and /2 are some functions of their arguments and their derivatives^. Essentially, we need a 
prescription for picking out the contributions that a particular hj makes to the linearized equations. 
As one might guess, such a unique decomposition is not always possible. 

The nature of any decoupling will of course depend on the symmetries of the background metric, 
which determines the functional form of C. Thus, one might also correctly guess that the spherical 
symmetry of the background Schwarzschild metric will allow for the separation of spherical harmonic 
modes, while the background Kerr metric, which obeys the far more restrictive axial symmetry, will 
not. As a preview of the next section we note that in the linear Schwarzschild analysis, three of the 
perturbed Ricci components take the form, 

\ CXD 

= ^f2{hi) Ui = , 

i=Q 

SO that a simple projection of {uj\ will yield f2{hj) — 0, which is one of the desired PDEs governing 
the metric perturbations. For the simpler case of decoupling modes then, the angular functions are 
known and the problem is then to determine the differential equations that the hi obey. Whereas, 
in the more general case the angular eigenfunctions must themselves be determined as part of the 
solution. For the case of the linearized Kerr metric, for instance, these angular functions can be 
shown to form a general Sturm-Liouville system, but they must be determined numerically, see e.g., 

i- 

Assuming that we have completely solved the linearized problem, i.e., once the angular eigen- 
functions and the PDEs relating the hi, are known, we can begin to quantitatively examine the 
second order perturbations. Consider truncating the main perturbation equation at 0{e^). It is not 
hard to convince oneself that the O(e^) terms occuring in the perturbed Einstein equations now 
yield the inhomogeneous system. 



where £ is exactly the same linear operator found in the first order analysis, and S is just the 
collection of quadratically occuring linear perturbations (assumed known from the linear analysis) 

^This is really just the statement that the linearly perturbed Einstein equations reduce to, R^^^ = 0, in vacuum. 
^Recall, for instance, that for Legendre functions, {ui\uj) = d6sin{6)uiUj oc Sij. 
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that are produced by keeping all of the O(e^) terms that appear in the perturbed Einstein equations. 
The physical significance of this inhomogeneous relationship, and the need for a clear understanding 
of the linear analysis, is now apparent. The quadratic, first order perturbation, terms act as an 
effective source term for the second order perturbations, allowing us to quantify the back reaction of 
the linear perturbations on the system. The determination of the PDEs for the second order pertur- 
bations now follows exactly the same procedure as in the linear analysis, but here the projections of 
the angular eigenfunctions will pick out components of the source terms from S. Note that without 
a clear analysis of the linear problem one would not know how to calculate what the projected 
components of the source terms are. While this may seem obvious, much of the linear Schwarzschild 
problem can be solved by simply assuming that the angular modes decouple, and setting various 
coefficients of cos"(0) equal to zero. If one chooses this, albeit simpler, solution technique for the 
linear problem the second order analysis is doomed. 

This completes our general discussion of MP methods, and we now move on to examine the 
particular case of perturbations to the Schwarzschild spacetime. The next section will demonstrate, 
with the aid of computer algebra, the decoupling of the linear Schwarzschild modes. As we shall see, 
the use of computer algebra allows the linear analysis to be examined in a thoroughly transparent 
manner. With the computer doing all of the work, we do not have to resort to invoking the standard 
function transformations, which mix the original perturbation functions, to simplify the calculations. 
The final result of the analysis presented in the next section will be a simple procedure for extracting 
the decoupled perturbation modes, of arbitrary order, for a Schwarzschild black hole. 

The Linear Schwarzschild Analysis 

The bulk of this section is the input / output from a GRTensorll session run on the Maple V platform. 
(The Maple worksheet itself will be made available for download from the GRTensor website.) While 
reading this section one should bear in mind what our objective is. Rather than actually deriving the 
standard linear Schwarzschild results, i.e. the Zerilli wave equation we simply want to show that 
the angular modes of the linearized Schwarzschild metric will decouple. In so doing we will develop 
a procedure to consistently extract a particular multipole from the perturbed Ricci tensor. Note 
that in what follows we consider only a single "z*'"' term from the general perturbation expansion 
presented in the last section, with its index suppressed. This is done for simplicity and clarity, and 
with no loss of generality due to the linearity of C. (The concerned reader can mentally place a 
summed index on all of the perturbation terms.) 

Our general approach to this problem begins by constructing the standard form of the covariant, 
linearly perturbed, Schwarzschild metric. 



For a discussion of the particular form of the perturbed metric that we use, which will be seen shortly, 
see e.g. [1]. The next step in this process is to calculate the exact contravariant metric tensor, i?"^, 
that is associated with the linearly perturbed covariant metric, gaf3, and then to linearize it w.r.t. 
e. The resulting object is the linearly perturbed contravariant metric tensor. 



(0) , (1) 



where, 



= d,^(^9 {9rr, 999, 9<p^, 9tt) = diag ((1 - 2M/r) 



-1 



,r2,r2sin2(0),-(l-2M/r)) . 




These linearized co/contra- variant forms of the metric tensor are the fundamental tools of the linear 
theory: With them one can calculate the linear perturbation to any tensor quantity by simply 
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following through the normal, unperturbed, calculation while dropping all of the non-linear e terms 
that appear. While this procedure can certainly produce calculations that would be intractable "by 
hand" , it is a very simple matter to instruct the Maple computer engine to follow this algorithmic 
procedure. To arrive at the linearized Einstein vacuum equations then, we simply instruct GRTensor 
to calculate the Ricci tensor from the linearly perturbed metric, while truncating the non-linear e 
ijGrms 

Of course, i?^"^ — 0, for the Schwarzschild spacetime, so that truncation of the non-linear e term 
yields the linearized Einstein equations, R^^^ = 0. 

To begin the Maple session we load the GRTensorll libraries, our new simplification/perturbation 
routinesn, and the linearly perturbed metric. 



> restart : 

> readlib(grii) : 

> grtensorO : 



GRTensorll Version 1.64 {R3) 
4 November 1997 
Developed by Peter Musgrave, Denis Pollney and Kayll Lake 
Copyright 1994 ~ 1997 by the authors. 



Latest version available from : http : // astro. queensu.ca/~ grtensor / 
To initiate help type 7 grtensor 



> mineO : 

> read 'myutils .mpl ' ; 

> qload(lpschw) : 
Calculating ds for Ipschw 



. . Done. (0.000000 sec.) 

Default spacetime — Ipschw 
For the Ipschw spacetime : 
Coordinates 
x( up ) 
^ [r9(j)t] 
Line element 



ds 



2 _ {l + eH2{r,t)u{e)) dr'^ 

1-2"^ 
r 



+ 2eHi{r,t)u{0) dr d t 



r'^ {l + eK{r,t)vi{9)) dO^ + sin( 6i )2 ( 1 + e K( r, t ) u( 6* ) ) d (/)^ 



-(^1-2— j (1 -eiJo(r,<)u(6')) df" 
Constraints = [e^ = 0, = 0, = 0, = 0, = 0, = 0, = 0, = 0, e^" = O] 



^ These routines are contained in the "myutils. mpl" file and are detailed in Appendix A. 
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These constraints will allow us to truncate the non-linear perturbation terms during calculation, 
which greatly reduces both the CPU and memory requirements. Note that the function u{d) ap- 
pearing in the metric represents an arbitrary Legendre function. (Although, it will actually remain 
an arbitrary function of 6 until we eliminate u" in terms of u' and u, using the Legendre equation, 
later in the calculation.) 

Our analysis begins with the calculation of the exact g"^ corresponding to the input metric, 
which we must then linearize in e. The linearization is accomplished using the linpert() routine, 
the details of which can be found in Appendix A. We are then left with, (7^/3 = g^^p + eg^^^, and. 
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9(0)"^ + ^9(1)°"^ 1 as the co/contra- variant components of the metric tensor^. 



> grcalc(g(up,up)) ; 

Calculating detg for Ipschw ... Done. (0.017000 sec.) 
Calculating g(up,up) for Ipschw ... Done. (0.067000 sec.) 

CPU Time = .067 



> gralter (g(up,up) .linpert , simplify .factor) ; 
Component simplification of a GRTensorll object: 

Applying routine linpert to object g(up,up) 
Applying routine simplify to object g(up,up) 
Applying routine factor to object g(up,up) 



CPU Time 



.133 



> grdisplay (g(up,up) ) ; 



For the Ipschw spacetime : 
Contravariant metric tensor 
g( up, up ) 
{-l + eH2{r,t)u{0)) (r-2m) 



,0,0,sHi{r,t)n{e) 



0^_^l±£KMW),o,0 



0,0, 



-l + eK{r,t)n{e) 
r2 ( cos( 6* ) - 1 ) ( cos( e) + l) 







eiri(r,t)u(0), 0, 0, 



r (l + £iJo(r,Ou(6')) 
r — 2m 



With the linearized co/contra- variant forms of the metric tensor in hand, it is a simple matter to 
calculate the perturbed Ricci tensor, or any other first order (tensor) quantities of interest for that 
matter. We now calculate the Ricci tensor, which amounts to calculating the perturbed Einstein 
tensor, applying the constraints as we go to kill off the higher order terms. 



> grcalcalter (R(dn,dn) , 13) ; 

Simplification will be applied during calculation. 

Applying routine Apply constraints repeatedly to object g(dn,dn,pdn) 
Calculating g(dn,dn,pdn) for Ipschw ... Done. (0.033000 sec.) 
Applying routine Apply constraints repeatedly to object Chr(dn,dn,dn) 
Calculating Chr(dn,dn,dn) for Ipschw ... Done. (0.050000 sec.) 
Applying routine Apply constraints repeatedly to object Chr(dn,dn,up) 

^Direct calculation shows that, g°'"'g~ifj = ~^ 0{€'^), as required for consistency at 0{e). 
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Calculating Chr(dn,dn,up) for Ipschw ... Done. (0.117000 sec.) 
Applying routine Apply constraints repeatedly to object R(dn,dn) 
Calculating R(dn,dn) for Ipschw ... Done. (0.733000 sec.) 



CPU Time = .933 



Up to this point the angular function u{9) has been completely arbitrary, except that one would like 
it to be a member of a complete set. We now fix u as a Legendre function by making the following 
substitution in Rap: u" — — cot{9)u' — n{n + l)u = 0. (For simplicity we set j — n{n + 1).) 



> grmap(R(dn,dn) , subs ,dif f (u(theta) , theta$2)= 

> -cos(theta) / sin(theta) *dif f (u(theta) ,theta) -j*u(theta) , 'x' ) ; 
Applying routine subs to R(dn,dn) 



We can now examine the seven non-trivial components of the perturbed Ricci tensor: ^rt\ 



, and the four diagonal components. The component is^ 



> bl : =hcollect (expand (coeff (collect (gr component (R(dn, dn) , [r ,theta] ) .epsilon) , epsilon, 1) ) , 

> ufuncs jlperts) ; 

[^2 {r^2m)r 2 {r-2m)r 2 \dr '7 2 \dr ^ ' ' 

'2 r-2m ) [de''^^^ 

(The hcollect() routine, described in Appendix A, is a hierarchical collection procedure, "ufuncs" 
is the set of u and its derivatives, and "Iperts" is the set of linear MP and their derivatives up to 
second order.) If we considered the full series perturbation mentioned above, we would have an 
infinite sum of these expressions, each with index i. Noting that the u'^{9) form an orthogonal set 
however, we could uniquely select out any desired term in the imagined summation by taking the 
appropriate projection on the linearized expression (and any potential source term). Examining i?!^^ 
we similarly find: 



> cl : =hcollect (coeff (collect (gr component (R(dn,dn) , [t , theta] ) , epsilon) , epsilon, 1) , 

> ufuncs , Iperts) ; 

Here again we simply need to project with u'^{9) in order to extract the "i*'"' term from any summed 
expression. We next find that R^W rI]^ , and iJ^j"* can all be expressed as: f{hi)ui{9) — 0. For 
these three components we need only project with the Ui function itself. 



> dl : =hcollect (expand (coeff (collect (gr component (R(dn, dn) , [t ,r] ) , epsilon) .epsilon, 1) ) , 

> ufuncs , Iperts) ; 

l2 r2 \drdt J r {r-2m)r / ^ ' 



^Although = 0, we will still take the time to collect only the linear e terms as a matter of good practice. 
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> el : =hcollect (expand (coeff (collect (gr component (R(dn,dn) , [r ,r] ) , epsilon) ,epsilon, 1) ) , 

> uf imc s , Ipert s ) ; 



\2r{r-2m)^2 (r-2m)2 r-2m ^ 2 \dr^ 

K{r,t)) I ^ ^^rHo{r,t)) m , 1 (-3m + 2r) {j: H,{r,t)) 



dr'^ ' / 2 r(r — 2to) 2 r {r — 2 m) 

(_3^ + 2r) (|:K(r,i)) m(fffi(r,t)) 



r(r — 2m) (r — 2m)^ 



> fl :=hcollect(expand(coeff (collect (grcomponent(R(dn,dn) , [t ,t] ) , epsilon) .epsilon, 1)) , 

> uf unc s , Ipert s ) ; 



1 {r-2m)^ (^^Ho{r,t)) i {r - 2 m,) H^,{r,t) j 

2 r2 "^2 r3 

1 (2r-m)(r-2m) {-§^Ho{r,t)) 1 m(r-2m) {^H^jr^t)) 

2 2 r3 
(r-2m)m(|:K(r,0) (-3m + 2r) gi(r,0) \ 

+ ;3 + ;:2 j^y^) 

So far all of the perturbed Ricci components have had a very simple angular dependence in terms of 
cither, u{6), or its derivative. The last two non-trivial components of i?^^^ do not, however, exhibit 
this simplicity. 

> al:=hcollect(coeff (collect (grcomponent(R(dn,dn) , [theta.theta] ) .epsilon) .epsilon, 1) . 

> ufuncs.lperts) ; 

al ..= [-{j + 2)H^{r,t) + -{-2+j)K{r,t)+ ^ 



2 ' -V ' / ■ 2 ' ^ ' ' ■ 2 r-2m 

-\ (|^K(r,0) {r-2m)r-\H,{r,t)j + \ (^^^H,{r,t)^ (r-2m) 

+ \ (^^2(r,0) (r-2m)-(-3m + 2r) (^i^K(r,o) " (^^i(^'O) <0) 

( 1 cos(g)go(r,t) 1 cos{6)H2{r,t) \ f d_ 
V 2 sin(e) ^2 sin(e) y'V^^'' 

Here we seem to run into trouble; our expression has both u'^ and Ui terms. The projection {ui\Uj) 
will in general be non-zero. Thus we could not decouple the modes if this expression were summed 
over. But when we combine this expression with i?^^^ we will find a, seemingly fortuitous, resolution. 

> gl : =kf actor (hcollect (coeff (collect (gr component (R(dn.dn) . [phi .phi] ) .epsilon) .epsilon. 1) 

> .ufuncs .Ipert s) . (cos(theta) -1) * (cos (theta)+l) ) ; 

gl := (cos(6»)-|-l)(cos(6l)-l) i- ^ (2H2{r,t)r-4:H2{r,t)m-2K{r,t)r + 4:K{r,t)m 
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+ K(r,Ojr-2K(r,Oim + r3 ( ^K(r,0 - _K(r,0 



+ 4 (|:^fo(r,o) (|:if2(r,0) ^2 - 4 (|:i^2(r,o) mr + 4 (|:^f2(r,0 ) m2 

+ 14 (|:K(-.0) "^--12 (|:K(-.0) m2-4 (|:K(r,i)) r2 - 2 ) r2 



+ 4 (^^ifi(r,i)) rm^u{e)/{r-2m) 

1 sin(g)cos(g) {H2{r,t) - Ho{r,t)) (^u(g))' 



2 (cos(6')-l)(cos(6') + l) / 

Here again we have the potential for mode couphng from the presence of both u'^ and Ui terms. If, 
however, we add/subtract these last two components a remarkable thing happens: 

> hcollect (al+gl/ sin(theta) ~2,uf uncs ,lperts) ; 

l-ij + 4)H,ir,t) + {-2+j)K{r,t)+ ^ _ K(r,t)J (r-2m)r 

-^Ho{r,t)j+(^-^Ho{r,t)^ {r - 2 m) + (J;^ H2{r,t)^ (r-2m) 
-2(-3m + 2r) (^|-K(r,o) -2 (^i?i(r,t)) rju(0) 

By adding these two expressions we have eliminated the u^ terms, and can therefore select any "i*'"' 
term by simply projecting with Ui. Similarly: 

> factor (hcollect (al-gl/ sin(theta) ~2 ,ufuncs ,lperts) ) ; 

2 sm{0) 

While this subtraction seems to again leave us with the potential for mode coupling, we notice that: 

i{i + + 2 cot(6')M-(6') = -(1 - x^)(fu^/dx^ , 

where x = cos{9). Recognizing this as the Pf{x) associated Legendre function, and recalling that 
these also form a complete orthogonal set, we can, therefore, uniquely select out any term in such a 
series by projecting with: 

i{i + l)ui + 2 cot{e)dui/d0 . 

The decoupling of the linear problem is thus reduced to a prescription for performing seven 
projections, 

from which one obtains seven PDEs that govern the MP functions. Here i = 1..7, Fi produces a 

linear combination of its arguments, and Gi is at most a function of the given argument and its 
first derivative. We can summarize the details of this decoupling procedure in the following table 
that relates the seven linear combinations of the perturbed Ricci components to their associated 
projection functions. 
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Perturbed Ricci components 


Z-mode Projection function 


6Rre , SRet 


dPi/de 


SRrr , SRrt , SRu 


Pi 


SReg + SR^^/am^{9) 


Pi 


SRgg - SR^^I sin"^(6') 


l{l + l)Pi + 2cot{0)dPi/d9 



Table 1. 



This table is the central result of this section, and is a kind of Rosetta stone for the Schwarzschild 
perturbation analysis. It allows one to consistently and algorithmically decouple the angular per- 
turbations, regardless of the order of the perturbation expansion^. Thus the SR^fs can be any order 
of the perturbed Ricci tensor, R^^^, or even the total perturbation, J2n ^"-^o/? ■ 

While we have managed to decouple the perturbation modes, we should note that there is a great 
redundancy in the system of seven PDEs that are generated by these projections. The remainder 
of the solution to the linear problem is, essentially, the elimination of this excess information. As a 
simple example of this redundancy, and as a prelude to the results of the next section, we note that 
the "96 — (fxj)" projection in the linear analysis gives, 

/ desm{9)P^{9) {SRee - SR^^/sm^{9)) = Ho{r,t) = H2{r,t) . 

At the linear level then, Hq and H2 are degenerate^ 

The reduction of the system is typically accomplished by making the Zerilli transformation, 

(2r^ — 6rm — 3m^) , , r'^r](r,t) dK(r,t) ^.ir^+rm + m^) , . , . 
[r — 2m)[2r + 3m) r ~ 2m ot r^(2r + 6m) 

With these substitutions the redundancy can now be expressed as: r]{r,t) = (1 — 2m/r)drX- Using 
this result to eliminate rj throughout the system one then finds a "wave equation" in the single new 
function x(r, i). This is the celebrated Zerilli "wave equation": 

-qT^ Q^ + yd-)xir,t)=0. 

Here, r* — r + 2Mlog(r/2M — 1), and the effective potential, V, has an explicit dependence on the 
angular mode number, I. Considering our discussion of perturbation theory from the last section, 
we expect that the second order Zerilli wave equation will be the same as that of the linear case, but 
with the addition of a new source term. Finding the form of this source term, which is not unique^, 
is the central problem of the second order calculation. Calculating the source term will, ultimately, 
be the goal of the next section. 

As the final comment of this section, we invite the interested reader to download the Maple 
worksheet that produced the linearized results presented here. One is then free to reproduce these 
results and experiment with the tools at hand. In particular one should try to replicate the linearized 
equations without the aid of the termsimp() and hcollect() routines. 



The Quadratic Schwarzschild Analysis 

We found in the last section that the non-trivial components of the linearly perturbed Ricci tensor 
grouped into five components, and two linear combinations of components, each having an asso- 
ciated projection function. Making use of this linear analysis we can now solve the second order 

® This is a direct result of the fact that higher order perturbations will simply produce the same seven equations, 
now in terms of the higher order functions, along with possible source terms. 

'^We will find that this degeneracy is broken at the second order level by the presence of "source terms" . 
®See, for example, [hi, for a discussion of the gauge dependence of the source term. 
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Schwarzschild problem. Here the second order perturbations appear exactly as the linear ones have, 
but with the addition of "source" terms that are quadratic in the (assumed known) linear perturba- 
tions. The key here is that while the quadratic source terms will in general be a combination of an 
arbitrary number of multipoles, we now have a unique prescription for determining their contribution 
to any desired second order multipole perturbation. 

From the point of view of the second order analysis, the form of the linear perturbation functions 
are, apart from consistency, simply a matter of choice. The recent work that has appeared in the 
literature assumes that the linear perturbations have only a quadrapole component, see, e.g. j^]. 
This restriction can be justified on physical ground by arguing that most of the spectral power will 
be concentrated in this mode and that the full source could be reconstructed by simply summing over 
the actual multipole modes of the linear problem. But perhaps the best rational for this restriction is 
that there exists a particular class of interesting problems, see, e.g. Q, for which the only non-zero 
multipole is, at the linear level, the quadrapole. For this case then, the treatment is exact, and 
it is this scenario that we will study. The remainder of this section will outline the calculation of 
the second order quadrapole results, i.e., g^^^, and, g^^^, will both be pure quadrapole moments. 
In solving this particular problem then, we are examining how the first order quadrapole moment 
seeds the second order quadrapole moment. We will return to discuss this point in our conclusion, 
but for now one should simply notice that the following analysis could easily be repeated with any 
two given multipoles. Further details of the calculation can be found in Appendix B, and as in the 
last section, the complete Maple worksheet that produced these results can be downloaded from the 
GRTensor website. 

The quadratically perturbed covariant metric for the "quadrapole-quadrapole" perturbation is a 
simple extension of the linear case. 



(^1 + i {eh2{r,t)+e^H2{r,t)) (3 cos( )2 - l)^ dr^ 



r 



+ {ehi{r,t)+e^Hi{r,t)) {3cos{df - l) dr d t 
+ (^^ + 1 (£k(r,0 +e^K(r,0) (3 cos( 6* )2 - l)^ d9 

+ r^sm{9)^ (^^ + \ {ek{r,t) + e^K{r,t)) (3 cos( 6* )2 - l) 

- (l-2^) (^1-i {ehoir,t)+s'Ho{r,t)) (3 cos( )2 - l)^ 



dt' 



Constraints = [e^ = 0,s^ = 0, = 0, = 0, = 0, = 0, e^" = 0, = O] 



As can be seen from the form of the metric, the lower case functions are the first order perturbations, 
and the upper case functions are the second order perturbations. 

The calculation now proceeds exactly as in the linear case, except that we now expand the 
contravariant metric to second order and truncate higher than O(e^) terms during the calculation 
of the Ricci tensor. One can gain an immediate appreciation for the increase in complexity of this 
problem over the linear case by examining the size of the second order Ricci tensor. 

R r r ~ 86717 words. Exceeds grOptionDisplay Limit 
RrO— 28579 words. Exceeds grOptionDisplayLimit 
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Rrt — 53004 words. Exceeds grOptionDisplay Limit 
R g = 87057 words. Exceeds grOptionDisplayLimit 
R t = 20438 words. Exceeds grOptionDisplayLimit 
R 4> (p = 72053 words. Exceeds grOptionDisplayLimit 
R 1 1 = 82672 words. Exceeds grOptionDisplayLimit 



As one might expect, the quadratically perturbed Ricci tensor is substantially larger than its linear 
counterpart, and it is here that our new simplification routines will become indispensable. Of course, 
what we have really calculated here is, 

Ra0 - Ra(3 + «^Q/3 + « ^ap ' 

where iJ^^"^ = 0. Given that the Ricci tensor now contains terms that can be either linear or quadratic 
in e, care must be taken when extracting a particular order of perturbation. 

In order to determine the PDEs that govern the second order perturbations we simply follow 
the procedure set out in the linear analysis. Here we select only the components of the Ricci 
tensor, and we perform the I = 2 projections. This projects the squared quadrapole of the linear 
perturbation onto the second order quadrapole]^ The second order perturbation functions will, 
necessarily, appear in these PDEs in exactly the same manner as the linear perturbation functions 
appear in the last section. The only new feature, which lies at the heart of the second order theory, 
is the manner in which the quadratically occuring linear perturbations appear. We will not show the 
raw expressions for most of these PDEs as many of their source terms are quite large. Instead, we 
give the detailed expressions for two of the simpler equations, both to illustrate how these equations 
can differ from their linear counterparts and to make contact with the results presented in |l|. The 
results of the "6*0 — 00" and projections yield, 

H2{r,t) = Ho{r,t) - l^hQ{r,tf + ]^hi{r,tf , 

and, 

„mi?i(r,t) (r-2m) (|:ii"i(r,t)) f d , A (3^^, A i fd, , 



respectively. The degeneracy between Hq and H2 is therefore broken at second order by the appear- 
ance of "source" terms. Of course, one can still eliminate H2 from the analysis using this first of 
these expressions. As with the linearized case, the second order perturbations form a very redundant 
system. The effective source terms, of course, inherit all of the redundancy of the linear problem as 
well. The completion of our second order calculation will be the elimination of this redundancy, d 
la Zerilli, in a manner that will allow us to reproduce the results of [Q. 

(2) 

The derivation of the first of the Zerilli equations proceeds by eliminating Hq from Rrr using 
(2) 

the i?^/ component of the Ricci tensor. One must then eliminate the redundant source terms, 
a procedure that is not unique. Here we choose to preferentially eliminate certain source terms to 
reproduce the results in ||l[ . The details of this process are not very enlightening and are relegated to 
the Maple worksheet that can be found at the GRTensor website. The result of these manipulations 
is essentially the raw result for the r]{r, t) expression which can be found in Appendix B. One should 
note that the quadratic nature of the linear perturbations that appear in these types of expressions 
make them almost impossible to simplify with the standard Maple routines. To overcome this 



^In the case of a more general problem, i.e., one that had more than a single second order perturbation, this 
procedure would also decouple the second order multipoles and accomplish the associated projections. 
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problem we use our new mcollect() routine (see Appendix A), after successive application of which, 
we find, 

(r-2m) (|:x(r,0) 1 / f d , . 

r]{r,t)^ - -{r -2m) \~ ho{r,t) I — /io(r,i) 

+ 2 (^/^o(r,t)) k(r,0+r (|:k(r,0) /.o( r, o) + 2 k( r, t ) (|k(r,0 

+ 4^^(-'^)'-(-'^)^2"-^^^-'^^(^"(-'^^)+MM) f^MM) 

+ 2 
+ 2 



dt 

{r-2m)h^{r,t)ho{r,t) (ik(r,t)) (f /ii(r,t)) 



;f k(r,0) feo(r,t)m 




r — 2 m 

This result, which appears here exactly as outputed from the Maple session, agrees with the result 
found in [1] . The new source terms that appear in this second order result add a considerable degree 
of complexity to the problem, as one must use this result to eliminate all of the r/ dependence, which 
includes up to third order derivatives, from the analysis. One should also note that without the 
use of our new simplification routines this analysis would not have been as seamless: Ideally, when 
performing computer algebra calculations, one should never have to resort to a "by hand" calculation 
as this increases the chance of error by a very large factor. Our new simplification routines were 
designed precisely to avoid the necessity of any "by hand" calculations. 

The development of the second order Zerilli "wave equation" follows a similar procedure to this 

f21 * ' (2) 

last result. We first take the time derivative of , eliminate Hq and Hq using -R^ , make the Zerilli 

(2) 

substitutions, and substitute for rj from our first Zerilli result. We then repeat this process for R^g 
and add the resulting equations such that the x' term vanishes. The result of these operations is 
essentially the raw form of the result we are seeking. Having 421 terms, this expression is, however, 
rather large. 

In order to make contact with the results presented in |^ , we now work with the "renormalized" 
perturbation function, defined by, 

2 r'^k(r,t) (^k(r,t)) 2 

x{r,t)=C{r,t) + - ^ ' ^ >' +-k{r,t)' 

7 2 r + 3 m 7 



This choice of transformation reduces the number of terms by over a factor of three. The next 
problem we must tackle is to eliminate all of the redundant information contained in the source terms 
for In principle this should be relatively simple, and certainly straight forward; in practice it is 
neither. The basic necessity is to use all of the PDEs found in the linear analysis to eliminate as many 
of the "higher order" derivatives as is possible, and to preferentially eliminate any hi dependence. 
The rational for this is that, to linear order, all of the radiation information is determined by the 
Zerilli function, tjj, which depends only on k and ho. One therefore expects that the source terms 
can be expressed solely in terms of these two functions, and ultimately, solely in terms of ip- Our 
approach is thus to express /iq in terms of V' and k once the hi dependence has been eliminated from 
the source terms. 

V'(r,t) - irk(r,<) j (2r + 3m) 



ho{r,t) = 3^ + jr-kir,t 

r [r — Zm ) \or 
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After performing all of these eliminationspl, making the above substitution for ho, and, further 
simplifying the result to eliminate the k dependence, we arrive at the unsimplified form of our 
final result, which is still rather large at 154 terms. To simplify this expression we once again use 
a repeated application of mcollect(). Here we must collect and simplify the 17 different types of 
quadratic V' terms that occur. Making one last set of substitutions, /i — (r — 2m), and, A = (2r+3m), 
and extracting the source terms from the resulting expression yields: 
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_ 12 



r4 ^3;^ 

( 2 + 4 m + 9 r TO^ 6 ) V'( r, t ) %2 



( 112 + 480 TO + 692 r"^ to^ + 762 r"^ to"^ + 441 r to* + 144 to* ) r, i ) \-^ip{r,t) 



d 



dt 



, a 
1 lat 



1 (18r3 -4r2TO-33rm2 -48to3) (^^-(^0) 
1 (12r3 + 36r2TO + 59rTO2 + 90to3) 



12 



( 2 + 9 TO + 6 TO^ — 2 to"^ — 15 r to* — 15 ) ^^{r^ty 



^ (r^+TOr + TO^) (f V;(r,t)) %1 ^ 

^2 



(32r* + 88r*TO + 296r3TO2 + 510r2TO3 + 561rm* + 270TO5)?/'(r,i) ^^^(^,0^ / 



1 ^('^^O) (o^V^('^^O) (2r^ -TO^) (^V,(r,0) %2 



3 fj, X 

{8r^ + 12mr + 7m^)^P{r,t)%l 1 (3r_7 to ) (|: V( r, i)) %1 

^ 3 



mVX^O 4 (3r2 + 5TOr + 6TO2) (|:^/>(r,0) %2 1^A%2^ 

^ 3 3 



1 A%r 



3 



/A 



%2 :=^V'(r,0 



This is our final result, and also appears exactly as outputed from the Maple worksheet. We have 
derived the effective Zerilli source term that dictates how the linear quadrapole perturbation seeds 
the second order quadrapole perturbation for a Schwarzschild black hole. This expression is almost 
identical to that presented in ||l[. The difference between the above expression and the previously 
published result is that the term that appears in j^] only has factor of l//i^, while we have found 

^"These details can be found in the worksheet at the GRTensor website. 
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that it has a l//i'^ dependence. 

To conclude this section we again invite the interested reader to download the full Maple work- 
sheet from the GRTensor website. One can then reproduce these calculations and, of course, exper- 
iment. One could, for instance, alter the worksheet to calculate the source term for any particular 
pair of linear-quadratic multipole perturbations. 



Conclusion and Discussion 

This article has presented a detailed account of our calculation of second order perturbations to a 
Schwarzschild black hole. Our methodology was chosen to illustrate both the utility, and necessity, of 
employing computer algebraic methods to examine these types of problems. We began our analysis 
with an explicit demonstration that the linear Schwarzschild perturbations decouple, a result that is 
well known but that is not always clearly presented. We then used these linear results to examine the 
second order pure quadrapole case. These second order results confirm, with one minor exception^, 
the earlier work by Gleiser, et. al, in [Q. 

As mentioned earlier, there are interesting problems for which the quadrapole perturbation pro- 
vides an exact description of the system to first order. Given this, one may then wish to consider how 
the linear quadrapole seeds an arbitrary second order multipole perturbation. This more general 
result has been presented by PuUin |^ . (This result contains a source term that remains inconsistent 
with our results.) With our general method it would take relatively little effort to examine this for 
any second order multipole^ or even to reproduce the general result from Our methods, of 
course, are not limited to examining perturbations that have only a linear quadrapole contribu- 
tion. Indeed, given the computational speed with which our simple "pure quadrapole" case was 
solvedPI, there is no reason to believe that more realistic perturbations could not be examined. 
One could, for example, study how the inclusion of the first few multipoles at the linear level affect 
any given second order multipole. One might even consider higher order calculations. Given the 
minimal computational requirements of our second order calculations it is likely that third order 
Schwarzschild perturbations would be manageable. While the utility of a third order calculation is 
certainly questionable, the potential for its investigation does exist. 

The largest obstacle to such calculations, and even to our present calculation, is finding the correct 
strategy for eliminating the redundancy in the linear perturbation equations. With some forethought, 
however, one could create a Maple routine that would examine a system of redundant PDEs, such 
as the ones we have encountered, and automatically generate a list of simplifying equations. One of 
these expressions might, for example, allow one to eliminate a second order derivative in terms of a 
linear combination of first order derivatives. Such a routine would certainly go a long way towards 
making these kinds of calculations much easier. 

With the methodology of Schwarzschild perturbations on such firm ground one is quite naturally 
led to carrying this analysis over to the calculation of perturbations to the Kerr metric. As mentioned 
in the introduction, this approach fails almost immediately when applied to the Kerr spacetime. The 
reason for this is simple: The background Kerr metric is not spherically symmetric, which results 
in perturbation equations that are not separable and that can therefore not be decoupled. One is 
free to try and find a coordinate system in which the equations become separable, but one should 
be forewarned that this is not a rewarding pursuit. None of the common coordinate systems of the 

^'^One, out of seventeen, of our Zerilli source terms differs from the previously published results. 

^^Gleiser, Q], for example, has considered the linear quadrapole in conjunction with the second order monopole 
perturbation. This is essentially an investigation of how the second order mass perturbation is affected by the first 
order gravitational radiation. 

^^The calculations presented in Appendix B were performed within a Maple worksheet running on a 600 MHz 
Alpha. The total CPU time was just under two minutes. 
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Kerr spacetime produce separable perturbation equations. Unless one can guess, or derive, such a 
coordinate system, this approach to the Kerr problem is a dead end. It thus appears that one must 
appeal to the tetrad solution of the linear Kerr problem for guidance. 

The central result of the tetrad analysis is the Teukolsky equation (see, e.g., Q), CT'<Pi = 0, where 
Ct is a complex linear operator and the V'i are the ^pQ and 1P4/ p'^ scalars from the NP formalism. All 
of the perturbation information is therefore encoded in just two complex scalars, so that one expects 
that there are only four independent MP functions. Our approach to this problem is to try and 
reverse engineer the NP solution and express the metric perturbations in terms of the perturbed NP 
quantities. If this can be done we will be assured that the resulting linear perturbation equations 
will decouple. By examining how the combinations of the separable eigenfunctions appear in the 
resulting linearized Ricci tensor one could then develop an algorithm, similar to our Schwarzschild 
result, to decouple the linear perturbations of the Kerr spacetime. If this linear problem can be 
solved with "reasonable" computing resources then one should be able to repeat the analysis for 
second order Kerr perturbations, just as we have done for the Schwarzschild case. 

Our preliminary calculations in this effort indicate that even the linear Kerr problem can gen- 
erate extremely large and complex expressions. Although much of this complexity disappears after 
extensive simplification. Despite the size of some of these intermediate objects, we are confident 
that the linear Kerr problem can be solved by this "reverse engineering" method. Our work is now 
proceeding in this direction, and we expect to have definite Kerr results shortly. Our findings here 
will be the subject of Paper II. If the linear Kerr calculation requires a large percentage of our full 
computing resources then it is likely that our findings for the second order Kerr case will simply be 
that a full solution to the problem will have to wait for the next generation of computer/computing 
engine. As our final comment we wish to note that the difficulty and complexity of the linear Kerr 
problem can be truly surprising. 
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Appendix A 

This appendix contains some of the details of our new Maple routines that were used in our calcu- 
lations. Most of these routines are fairly small, and did not require a significant amount of coding. 
(There are, however, some notable exceptions.) Our exampleg|^ will make use of the following two 
polynomials that are similar to those found in our calculations: 

> pi : = (r ~2+M*r+M~2) ~2* (r-2*M) -2/r ~2/ (2*r+3*M) ~2; 

{r^ + M r + M^y {r - 2 Mf 



pi :-- 



r^{2r + 3My 



> p2 : = (2*r"2-6*M*r-3*M"2) '2* (r-2*M) '2/r/ (r-2*M) "2/ (2*r+3*M) ; 



p2 :-- 



r{2r + 3M) 



• termsimp(ea;pr): The termsimp() routine applies factor (simplify ()) to each term in a type + 
or * expression individually, and then reconstructs the final expression. This can be very useful if one 
needs to simplify an expression after having applied the collect() routine, as Maple's simplify() 
and collect() tend to be inverse operations. 

e.g. 

> temp : =expand(pl) *H [0] +expand(p2) *H [1] ; 

temp := — 2 — — ,„ — 2 ■ 



(2r + 3M)2 (2r + 3M)2 (2r + 3Af)2 (2r + 3M)2 (2r + 3M)2 

r(2r + 3M)2 r^{2r + 3M)^J 

r^M rM"^ \ „ 

24 TTTT + 24 -— + 36 — + 9 — — — Hi 



2r + 3M 2r + SM 2r + 3M 2r + 3M r(2r + 3M) 

> termsimpCtemp) ; 

(r - 2 M )^ {r^ + M r + M^f Ho {2r^ - Q M r - i M"^ f Hi 
r2 (2r + 3M)2 r ( 2 r + 3 M ) 



We can now compare this with the result of applying the Maple simplify() routine 
> simplify (temp) ; 

[Ho -2Ho M - H^ - 2 + 5 Hq r"^ M'^ + i Hq r + iHg + 8H 

- 36 Hi M - 24 Hi + 144 Hi + 126 Hi M* + 27Hir M' 

{2r + iMf 



This rather cumbersome expression is a result of the fact that simplify () does not "know" that the 
i/i's should be treated as the primary objects. The effect of using collect () (which would produce 
an expression similar to temp) and then termsimp() essentially tells Maple that the "collected" 
objects are the primary ones. 

^'^ While these examples may seem somewhat artificial to the uninitiated, one should bear in mind that the result 
of a large algebraic calculation is typically an almost fully expanded expression that must be simplified. 
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• hcollect ( expr, list_one, list_two): The hcollect() routine is a hierarelial eolleetion routine. Tlie 
collect() routine is first applied to the expression with list_one as an argument. The resulting 
expression, which has been factored w.r.t. the elements of list_one, is then collected w.r.t. the 
elements of list-two. The utility of this function is that it is designed not to undo any of the first 
collection, which is the typical result of applying successive collect() calls with different lists. This 
routine automatically calls termsimp() when reconstructing the final expression. 

e.g. 

> temp : =expaiid(pl* (u(r) + (r-2*M) /r*dif f (u(r) ,r) ) *H [0] + 

> p2*(u(r)+(2*r+3*M)*diff (u(r) ,r))*H[l]) ; 

^'^o'^(0 , ^'Ho{§^n{r)) Hp n{r)) M ^ r'HoMnjr) 

^' (2r + 3M)2^ (2r + 3M)2 (2r + 3M)2 (2r + 3M)2 

r2goM2 (|:u(r)) r^HoM^ujr) rHpM^n{r) H^M^ (|:u(r)) 
(2r + 3M)2 (2r + 3M)2 (2r + 3Af)2 ^ (2r + 3M)2 

HoM^ujr) HqM^ (|:u(r)) goM-^u(r) Hp (|:u(r)) 

(2r + 3M)2 r(2r + 3M)2 ^ r{2r + 3Mf r2(2r + 3M)2 

HoM^ujr) HqM' {-§,u{r)) H,u{r) H, (^u(r)) 

r2(2r + 3M)2 r3(2r + 3M)2 2r + 3M 2r + 3M 



r 



^H,{Iu{t))M ^^r'H.Mujr) ,,r^H,M' (|:u(r)) , ^^rgrM2u(r) 



_ 3g ^ \ar V ^/ 24 ^ ' — 24 ^"'^ ' + 24 

2r + 3M 2r + 3M 2r + 3M 2r + 3M 

r//ii\/'' (-f u( /■)) iJii\/'u(r) (^u(r)) IhM^Axir) 

+ 144 1 V,;r V _^3g_L L^_^i26^ " ) +9 



+ 27 



2'r + 3Af 2r + 3M 2r + 3M r(2r + 3M) 

H^M^ (^u(r)) 



r (2r + 3M) 
> hcollect(temp,{u(r) ,diff (u(r) ,r)}-,{H[0] ,H[1]}) ; 

(r2 + Mr + M2)2(r-2M)3ifo {^2r^ -%Mr -"iW^f ( d 



r3(2r + 3M)2 r ) \dr 

/'{r-2Mf{r^+Mr + M'^fHQ {2r^ -6Mr -3M'^f HA 

r2(2r + 3M)2 + r(2r + 3M) j ''^''^ 

As with any "free code" , one should always perform some random consistency checks: 

> simplify(hcollect(temp,-[u(r) ,diff (u(r) ,r)>,{H[0] ,H[l]})-temp) ; 





• mcollect(ea;pr, {arg_one [, arg_two ]}): The mcollect() routine will collect the terms of an 
expression w.r.t. the product {arg-one argJ,wo). If argJ,wo is not supplied the collection is done 
w.r.t. {arg-one)"^. (This routine required a relatively large amount of coding. In particular, one has 
to very careful, from a coding standpoint, when the arguments can be derivatives of a function.) 
This routine automatically calls termsimp() when reconstructing the final expression. 

e.g. 

> temp:=expand(pl*(u(r)+(r-2*M)/r*diff (u(r) ,r))*(u(r)+(r-2*M)/r*diff (u(r) ,r))) ; 
^' (2r + 3M)2 (2r + 3M)2 (2r + 3M)2 (2r + 3M)2 
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(2r + 3M)2 (2r + 3M)2 (2r + 3M)2 r(2r + 3M)2 

{2r + ZMy r2(2r + 3M)2"^ r{2r + ZMy r4(2r + 3M)2 
r-4u(r)2 ^r4(|:u(r))\^ M6u(r)2 ^ r3u(r) (|:u(r)) M 



(2r + 3M)2 (2r + 3Af)2 r2(2r + 3M)2 (2r + 3M)2 

r3(|:u(r))^M r2M2u(r) (i|u(r)) u( r ) u( r )) 

(2r + 3M)2 ^ (2r + 3Af)2 (2r + 3Af)2 

M^u(r) (^u(r)) MSu(r) (^u(r)) M^u(r) (|:u(r)) 
r(2r + 3M)2 r2(2r + 3M)2 r3(2r + 3M)2 

> mtemp : =incollect (mcollect (mcollect ( 

> mcollect (temp, {diff(u(r),r)}) ,-Cdiff (u(r) ,r) ,u(r)}-) ,{u(r)}-) ,{diff (u(r) ,r)}) ; 

(r2+Mr + M2)2(r-2M)2u(r)2 {r"^ + M r + M'^f {r - 2 Mf Vi{r) (^u(r)) 
mtemp := h 2 ■ 



+ 



r2(2r + 3M)2 r3(2r + 3M)2 

(r2 + Mr + M2)2(r-2M)4 (|:u(r))' 



r4(2r + 3M)2 
And just to be sure: 

> simplify (mtemp-temp) ; 



• kfactor(ea;pr, fctr): The kfactor() routine simply forces Maple to pull the factor fctr out of the 
expression. The chief utility of this is largely cosmetic, although it can be useful when comparing 
expressions to published results. 

e.g. 

> temp:=pl*H[0]+p2*H[l] ; 

_ [r'^ + M r + M'^f (r - 2Mf Hq {2r'^ -QMr -iW^f 
^^"^^ r2(2r + 3M)2 + r(2r + 3M) 

> kf actor (temp, (r"2+M*r+M*2) "2* (r-2*M) *2) ; 
{r^ + Mr + M^f{r-2M) 



,2 , , ,^2^2. o,^^2^ Ho , (2r2-6Mr-3M2)2Hi 



r2(2r + 3M)2 r {2r + 3 M ) {r'^ + M r + M"^ )^ {r - 2 M y 



• linpert(ea;pr): The linpert() routine returns the first order Taylor expansion in e, about e = 0, 
of the supplied expression. While this routine could be made far more elaborate in terms of options 
and parameters, it was specifically designed as a single argument function so that it could be used 
within GRTensorll in a seamless manner. 

e.g. 

> temp:=l/sqrt(l+epsilon*f (r,tlieta,pM,t)) ; 

1 

temp := 



Vl + ef(r,0,(/.,O 
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> linpert (temp) ; 

l-ief(r,0,0,O 



• quadpert(e.Tpr): The quadpert() routine is similar to the linpert() routine, except that it 
returns the second order Taylor expansion of the expression in e. 

e.g. 

> quadpert (temp) ; 



